Анализ сортов сои

Правки

1. Введение + EDA

• перевели забытые переменные в факторные
• проанализировали в каких переменных встречаются пропущенные значения
• при поиске выбросов использовали пакет outliers, предположили, что выбросы это значения, которые отклоняются более чем на 1.5 межквартильного размаха от первого и третьего квартилей.

2. ANOVA:

• добавили осмысленности гипотезам
• провели оработку проблемы допущения с равенством дисперсий
• поправили графики
• провели численные тесты для проверки допущений

3. Линейные модели
4. Общие моменты

• добавили описание новых пакетов
• добавили меню для переключения между разделами
• подправили структуру
• обновили выводы

Используемые библиотеки

Новые библиотеки:

• plotly: визуализация данных (построение графиков)

• patchwork: объединение графиков и более гибкая настройки визуализации результатов ggplot2

• effects, sjPlot: проверка условий применимости смешанных линейных моделей

• purrr: упрощает сложные манипулцяии с данными, с помощью функции map в пакете возможно значительно сократить описание функции и переменных

Введение

Данная работа посвящена анализу сортов льна, выращиваемым в течение 5 лет на 2 разных локациях: на Кубани и в Липецке. Начнём анализ с изучения переменных, по которым была собрана информация:

## 'data.frame':    2943 obs. of  16 variables:
##  $ X                : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ id               : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ leaf_shape       : chr  "round" "round" "round" "round" ...
##  $ maturation_group : int  5 6 5 3 5 3 5 5 6 5 ...
##  $ lodging_type     : chr  "yes" "yes" "no" "no" ...
##  $ growth_type      : chr  "indeterminant" "indeterminant" "indeterminant" "semi_determinant" ...
##  $ flowering_group  : num  5 5 5 4 5 2 4 5 5 5 ...
##  $ pubescence_colour: chr  "light_gray" "light_gray" "light_gray" "light_gray" ...
##  $ corolla_colour   : chr  "white" "purple" "purple" "purple" ...
##  $ origin           : chr  "Uzbekistan" "India" "China" "China" ...
##  $ productivity     : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ vegetation_period: int  NA NA NA NA NA NA NA NA NA NA ...
##  $ protein_content  : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ oil_content      : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ site             : chr  "kub" "kub" "kub" "kub" ...
##  $ year             : int  2017 2017 2017 2017 2017 2017 2017 2017 2017 2017 ...

Заметим, что в исходном датасете есть пропущенные значения и категориальные переменные, которые для дальнейшей работы необходимо будет перевести в факторные.

Переменная	Описание	Вариации
id	Номер сорта	300 сортов
leaf_shape	Форма листа	lanceolate/round
maturation_group	Группа созревания	от 1 (раннеспелый) до 6 (позднеспелый)
lodging_type	Тип полегания	no/leaning/yes
growth_type	Тип развития	determinant/indeterminant/semi_determinant
flowering_group	Группа цветения	от 1.00 (рано) до 5 (поздно), шаг 0.5
pubescence_colour	Цвет опушения	gray/light_gray/light_tawny/tawny
corolla_colour	Цвет венчика	purple/white
origin	Происхождение	Страна
productivity	Продуктивность	г/м2
vegetation_period	Период от всхода до сбора урожая	дни
protein_content	Содержание белка в семенах	1-100%
oil_content	Содержание масла в семенах	1-100%
site	Место выращивания	kub (Кубань), lip (Липецк)
year	Год сбора	2017-2021

В выборке встречаются сорта родом из 31 страны.

##  [1] "Austria"     "Belarus"     "China"       "Czech"       "Estonia"    
##  [6] "France"      "Georgia"     "Germany"     "Hungary"     "India"      
## [11] "Italy"       "Japan"       "Kanaues"     "Kazakhstan"  "Latvian"    
## [16] "Lithuania"   "Moldova"     "Netherlands" "NorthKorea"  "Poland"     
## [21] "Portugal"    "Romania"     "Russia"      "Serbia"      "Slovakia"   
## [26] "South_Korea" "Sweden"      "Ukraine"     "USA"         "Uzbekistan" 
## [31] "Yugoslavia"

Переменные leaf_shape, lodging_type, growth_type, pubescence_colour, corolla_colour, origin, site переведём в факторные.

Maturation_group и flowering_group — в упорядоченные факторные.

Productivity, vegetation_period, protein_content, oil_content - оставим в первозданном виде.

data <- data %>%
   mutate(across(c(
    leaf_shape, maturation_group, lodging_type, growth_type, flowering_group, 
    pubescence_colour, corolla_colour, site, year, origin), as.factor)) %>%
  
  mutate(
    maturation_group = factor(maturation_group, 
                              levels = c(1, 2, 3, 4, 5, 6), 
                              ordered = TRUE),
    flowering_group = factor(flowering_group, 
                             levels = c(1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5), 
                             ordered = TRUE)
  )

Разведочный анализ данных

Рассмотрим данные с помощью функции summary

##        X                id             leaf_shape   maturation_group
##  Min.   :   1.0   Min.   :  1.0   lanceolate: 225   1: 117          
##  1st Qu.: 736.5   1st Qu.: 84.0   round     :2718   2: 468          
##  Median :1472.0   Median :167.0                     3:1017          
##  Mean   :1472.0   Mean   :166.4                     4: 828          
##  3rd Qu.:2207.5   3rd Qu.:249.0                     5: 378          
##  Max.   :2943.0   Max.   :330.0                     6: 135          
##                                                                     
##   lodging_type            growth_type   flowering_group   pubescence_colour
##  leaning: 162   determinant     : 936   3      :1449    gray       :  90   
##  no     :2655   indeterminant   : 729   5      : 414    light_gray :1701   
##  yes    : 126   semi_determinant:1278   4      : 387    light_tawny:1143   
##                                         3.5    : 234    tawny      :   9   
##                                         4.5    : 198                       
##                                         2      : 162                       
##                                         (Other):  99                       
##  corolla_colour     origin      productivity   vegetation_period
##  purple:2331    Russia :1557   Min.   :  9.9   Min.   : 72.0    
##  white : 612    Ukraine: 189   1st Qu.:128.0   1st Qu.: 96.0    
##                 France : 153   Median :173.0   Median :104.0    
##                 Kanaues: 126   Mean   :170.2   Mean   :103.9    
##                 China  : 117   3rd Qu.:210.0   3rd Qu.:113.0    
##                 (Other): 693   Max.   :362.0   Max.   :137.0    
##                 NA's   : 108   NA's   :1543    NA's   :1533     
##  protein_content  oil_content     site        year    
##  Min.   :33.00   Min.   :14.30   kub:1635   2017:327  
##  1st Qu.:39.50   1st Qu.:20.00   lip:1308   2018:654  
##  Median :41.20   Median :21.20              2019:654  
##  Mean   :41.07   Mean   :21.05              2020:654  
##  3rd Qu.:42.60   3rd Qu.:22.30              2021:654  
##  Max.   :51.50   Max.   :25.10                        
##  NA's   :1686    NA's   :1686

Опечатки

Поскольку в датасете довольно много категориальных описательных переменных мы проверили их на наличие опечаток и ничего не обнаружили.

## [1] lanceolate round     
## Levels: lanceolate round

## [1] determinant      indeterminant    semi_determinant
## Levels: determinant indeterminant semi_determinant

## [1] leaning no      yes    
## Levels: leaning no yes

## [1] gray        light_gray  light_tawny tawny      
## Levels: gray light_gray light_tawny tawny

## [1] purple white 
## Levels: purple white

## [1] kub lip
## Levels: kub lip

Пропущенные значения

Всего в датасете у нас 6556 пропущенных значений. На первый взгляд их довольно много. Подробнее изучим в каких столбцах они встречаются:

• origin содержит 108 пропущенных значений.
• productivity содержит 1543 пропущенных значений.
• vegetation_period содержит 1533 пропущенных значений.
• protein_content содержит 1686 пропущенных значений.
• oil_content содержит 1686 пропущенных значений.

При удалении всех пропущенных значений объем данных для дальнейшего анализа уменьшается в два раза, но при этом еще и полностью пропадают данные о познеспелых сортах (6), поскольку по ним нет измерений по количественным признакам. Кроме того, мы заметили, что за 2017 год в Липецке не проводилось измерений:

table(data$year, data$site)

##       
##        kub lip
##   2017 327   0
##   2018 327 327
##   2019 327 327
##   2020 327 327
##   2021 327 327

Мы можем позволить себе удалить все пропущенные значения, поскольку данных будет всё ещё достаточно для проведения статистического анализа:

data_without_na <- data %>%
  filter(if_all(everything(), ~ !is.na(.) & . != "" & trimws(.) != ""))

summary(data_without_na)

##        X                id             leaf_shape   maturation_group
##  Min.   : 141.0   Min.   :  1.0   lanceolate: 105   1: 54           
##  1st Qu.: 467.5   1st Qu.:117.0   round     :1090   2:260           
##  Median : 793.0   Median :179.0                     3:510           
##  Mean   :1078.5   Mean   :177.0                     4:259           
##  3rd Qu.:1705.5   3rd Qu.:239.5                     5:112           
##  Max.   :2826.0   Max.   :330.0                     6:  0           
##                                                                     
##   lodging_type            growth_type  flowering_group   pubescence_colour
##  leaning:  69   determinant     :346   3      :693     gray       : 38    
##  no     :1098   indeterminant   :282   4      :129     light_gray :685    
##  yes    :  28   semi_determinant:567   2      : 99     light_tawny:470    
##                                        3.5    : 91     tawny      :  2    
##                                        5      : 76                        
##                                        4.5    : 59                        
##                                        (Other): 48                        
##  corolla_colour     origin     productivity   vegetation_period protein_content
##  purple:942     Russia :687   Min.   : 12.0   Min.   : 72       Min.   :33.00  
##  white :253     France : 85   1st Qu.:132.0   1st Qu.: 96       1st Qu.:39.50  
##                 Ukraine: 81   Median :177.0   Median :104       Median :41.20  
##                 Kanaues: 80   Mean   :173.7   Mean   :104       Mean   :41.02  
##                 Austria: 45   3rd Qu.:213.0   3rd Qu.:113       3rd Qu.:42.60  
##                 Serbia : 36   Max.   :362.0   Max.   :137       Max.   :51.50  
##                 (Other):181                                                    
##   oil_content     site       year    
##  Min.   :14.30   kub:844   2017:173  
##  1st Qu.:20.00   lip:351   2018:401  
##  Median :21.20             2019:351  
##  Mean   :21.06             2020:161  
##  3rd Qu.:22.40             2021:109  
##  Max.   :25.10                       
## 

Выбросы

Проанализируем выбросы в численных переменных productivity, protein_content, oil_content. Стоит учесть, что это непростая задача, поскольку в выброке 300 разных сортов, обладающих уникальными признаками, потому выборка гетерогенна и сорта высоко/низкопродуктивные могут быть приняты как выброс.Поэтому данные для удаления стоит анализировать в совокупности с остальными признаками.

Будем искать с помощью пакета outliers. Предположим, что выбросы это значения, которые отклоняются более чем на 1.5 межквартильного размаха от первого и третьего квартилей.

1. Продуктивность

outliers <- boxplot.stats(data_without_na$productivity)$out

ggplot(data_without_na, aes(x = maturation_group, y = productivity, group = maturation_group)) +
  geom_boxplot(
    outlier.shape = NA, 
    color = "black", 
    lwd = 1.0
  ) +
  geom_jitter(
    aes(color = ifelse(productivity %in% outliers, "outlier", "normal")),
    width = 0.2, 
    height = 0, 
    size = 1
  ) +  
  scale_color_manual(values = c("normal" = "blue", "outlier" = "#c45824")) +  
  theme_minimal(base_size = 14) +
  labs(
    x = "Группа созревания",
    y = "Продуктивность",
    color = "Тип данных"
  ) 

2. Содержание белка и содержание масла

Можно заметить, что типы созревания, в которых встречаются выбросы по содержанию белка и масла совпадают (в групах 2,3,4 и 5). Кроме того, можно заметить, что выбросы, которые детектируются по большей части обратно-пропорциональны друг другу. Логично, что при высоком содержании масла будет низкое содержание белка. Поэтому мы решили не выбрасывать эти значения.

outliers_p <- boxplot.stats(data_without_na$protein_content)$out

plot1 <- ggplot(data_without_na, aes(x = maturation_group, y = protein_content, group = maturation_group)) +
  geom_boxplot(
    outlier.shape = NA, 
    color = "black", 
    lwd = 1.0
  ) +
  geom_jitter(
    aes(color = ifelse(protein_content %in% outliers_p, "outlier", "normal")),
    width = 0.2, 
    height = 0, 
    size = 1
  ) +  
  scale_color_manual(values = c("normal" = "blue", "outlier" = "#c45824")) +  
  theme_minimal(base_size = 12) + 
  labs(
    x = "Группа созревания",
    y = "Содержание белка %",
  ) +
  theme(
    legend.position = "none",  
  )

outliers_o <- boxplot.stats(data_without_na$oil_content)$out

plot2 <- ggplot(data_without_na, aes(x = maturation_group, y = oil_content, group = maturation_group)) +
  geom_boxplot(
    outlier.shape = NA, 
    color = "black", 
    lwd = 1.0
  ) +
  geom_jitter(
    aes(color = ifelse(oil_content %in% outliers_o, "outlier", "normal")),
    width = 0.2, 
    height = 0, 
    size = 1
  ) +  
  scale_color_manual(values = c("normal" = "blue", "outlier" = "#c45824")) +  
  theme_minimal(base_size = 12) + 
  labs(
    x = "Группа созревания",
    y = "Содержание масла %",
  ) +
  theme(
    legend.position = "none",  
  )

combined_plot <- (plot1 | plot2) +
  plot_annotation(
    title = "Поиск выбросов",
    theme = theme(
      plot.title = element_text(size = 14, face = "bold", hjust = 0.5)
    )
  )

combined_plot

3. Более изощренный поиск выбросов с помощью расстояния махаланобиса. Давайте попробуем проанализировать все численные данные по годам. Предположим, что порог для выбросов - 90%. Давайте посмотрим как распределяются все наши численные данные внутри каждого года. Будем смотреть в зависимости от года, поскольку в зависимости от погодных условий в течение года могут меняться и физиологические параметры растений.

# Визуализация многомерных выбросов по годам
library(ggplot2)

ggplot(numeric_data, aes(x = year, y = mahalanobis_dist, color = outlier)) +
  geom_boxplot(outlier.shape = NA, fill = "grey", alpha = 0.5) +
  geom_jitter(width = 0.2, height = 0, aes(color = outlier), size = 2) +
  scale_color_manual(values = c("normal" = "blue", "outlier" = "#c45824")) +
  theme_minimal(base_size = 10) +
  labs(
    title = "Многомерные выбросы по годам",
    x = "Год",
    y = "Mahalanobis Distance",
    color = "Тип выброса"
  ) +
  theme(legend.position = "right")

И несмотря на такой изощренный способ поиска, мы настаиваем, что лучше оставить значения, которые сильно отклоняются, потому что эти данные могут говорить о каких-то довольно ценных сортах.

Проверим распределения численных величин. С помощью функции map в пакете purr возможно значительно сократить описание функции и переменных, применив map(var, func)

## Warning: Removed 1543 rows containing non-finite outside the scale range
## (`stat_density()`).

## Warning: Removed 1686 rows containing non-finite outside the scale range
## (`stat_density()`).
## Removed 1686 rows containing non-finite outside the scale range
## (`stat_density()`).

В целом можно отметить, что численные величины распределены относительно равномерно (имееются несклько пиков), но не наблюдаем значительных отклонений.

Проверка гипотез

Первичный анализ корреляций

Знания биологии позволяют нам сказать, что содержание масла прямопропорционально связано с продуктивностью, тогда как содержание белка - обратно. Давайте посмотрим, насколько это воспроизводится на наших данных.

plot_1 <- ggplot(data_without_na, aes(x=productivity, y=oil_content))+
  geom_point(shape = 21, size = 1, fill = "#c45824", color = "black") +
  theme_minimal() +
  labs(
    x = "Продуктивность (г/м²)",
    y = "Содержание масла (%)"
  )

plot_2 <- ggplot(data_without_na, aes(x=productivity, y=protein_content))+
  geom_point(shape = 21, size = 1, fill = "#c45824", color = "black") +
  theme_minimal() +
  labs(
    x = "Продуктивность (г/м²)",
    y = "Содержание белка (%)"
  ) 

combined_plot <- (plot_1 | plot_2) +  plot_annotation(
    title = "Распределение данных по содержанию масла и белка",
    theme = theme(
      plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
    )
  )

combined_plot

Проверим статистичекую значимость кореляций продуктивности с различными параметрами:

Сначала посмотрим на нормальность наших выборок с помощью shapiro.test

Нулевая гиполеза H0 - выборки распределены нормально

Альтернативная H1 - выборки распределены ненормально

shapiro_var1 <- shapiro.test(data$productivity)
print(paste("p-value для productivity:", shapiro_var1$p.value))

## [1] "p-value для productivity: 0.00644333565970214"

shapiro_var2 <- shapiro.test(data$oil_content)
print(paste("p-value для oil_content:", shapiro_var2$p.value))

## [1] "p-value для oil_content: 1.36929282750539e-11"

shapiro_var2 <- shapiro.test(data$protein_content)
print(paste("p-value для protein_content:", shapiro_var2$p.value))

## [1] "p-value для protein_content: 1.38191116472042e-06"

Видим, что выборки распределены ненормально (p-value<0.5), применим корреляцию по Спирману (продуктивность и содержание масла):

cor_test_result <- cor.test(data$productivity, data$oil_content, method = "spearman")

print(paste("Коэффициент корреляции Спирмена:", cor_test_result$estimate))

## [1] "Коэффициент корреляции Спирмена: 0.514123548211428"

print(paste("p-value:", cor_test_result$p.value))

## [1] "p-value: 1.53882917165075e-85"

Продуктивность и содержание белков:

cor_test_result <- cor.test(data$productivity, data$protein_content, method = "spearman")

print(paste("Коэффициент корреляции Спирмена:", cor_test_result$estimate))

## [1] "Коэффициент корреляции Спирмена: -0.382382602249352"

print(paste("p-value:", cor_test_result$p.value))

## [1] "p-value: 6.19818358642783e-45"

Согласно шкале Чедокка, при абсолютном значении коэффициента корреляции Спирмена < 0.03 корреляция слабая, тогда как при от 0.5 до 0.7 заметная. В нашем случае мы обнаружили отрицательную корреляцию между продуктивностью и содержанием белков, тогда как продуктивность и содержание масел коррелируют более значимо.

Взглянем на коэфициенты, используя для этого пакет psych.

Прежде отберем в отдельный датафрейм интересующие нас количественные переменные

data_num <- subset(data, select = c(productivity, oil_content, protein_content))

cor(data_num, use = "complete.obs")

##                 productivity oil_content protein_content
## productivity       1.0000000   0.5475153      -0.4543480
## oil_content        0.5475153   1.0000000      -0.5929773
## protein_content   -0.4543480  -0.5929773       1.0000000

Наблюдаем отицательную корреляцию между продуктивностью, группой цветения и содержанием белков. Корреляция положительна между продуктивностью и содержанием масла. Всё как по учебнику.

Анализ различий между группами

Влияние места происхождения на продуктивность

H0 - место происхождения не влияет на продуктивность

H1 - место происхождения влияет на продуктивность

Посмотрим на боксплоты продуктивности в зависимости от места происхождения:

По построенным боксплотам видно, что среднее групп различается, поэтому гипотезу стоит проверить

model2 <- aov(productivity ~ origin, data = data_without_na)
summary(model2)

##               Df  Sum Sq Mean Sq F value Pr(>F)    
## origin        27  792464   29351   8.771 <2e-16 ***
## Residuals   1167 3905291    3346                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

По полученным графикам построенной модели, определим соблюдение допущений:

Остатки распределены нормально и линейно - допущение выполняется, однако есть проблемы с графиком дисперсий, попробуем численные тесты для проверки:

performance::check_homogeneity(model2)

## Warning: Variances differ between groups (Bartlett Test, p = 0.000).

p-value < 0.05, значит, что дисперсии между группами не равны.

Попробуем провести тест анова Уэлча, которая является модификацией теста классической анова не требующей равенства дисперсий:

model2 <- oneway.test(productivity ~ origin, data = data_without_na, var.equal = FALSE)
model2

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  productivity and origin
## F = 85.125, num df = 27.000, denom df = 31.038, p-value < 2.2e-16

p-value - значительно ниже 0.05, поэтому мы можем отклонить нулевую гипотезу и предположить, что место происхождения влияет на продуктивность сортов.

Однако теперь невозможно использовать классические post hoc тесты, потому что они необходимы для моделей с равенством дисперсий. Все-таки попробуем создать модель несмотря на разницу дисперсий между группами и уже потом провести poc-hoc тесты.

model_disp <- aov(productivity ~ origin, data = data_without_na)

Теперь попробуем выяснить, между каким группами возникают отличия. Для этого проведем тест Тьюки (Tukey HSD), который более строгий и выделяет следующие отличия между странами:

Тест Тьюки более строгий и выделяет следующие отличия между странами:

• China - Austria: p = 0.0009

• Japan - Austria: p = 0.0010

• USA - Austria: p = 0.0092

• China - Belarus: p = 0.0013

• Japan - Belarus: p = 0.0009

• Russia - China: p = 0.0000

• Ukraine - China: p = 0.0035

• Czech - China: p = 0.0110

• Japan - France: p = 0.0216

• USA - Czech: p = 0.0223

Отметим, что отличия между местами произрастания замечены в странах находящихся в разных частях света и соотвественно различных климатических зонах: можно выделить страны СНГ (Россия, Беларусь, Украина), страны Азии (Китай, Япония), страны европейской части (Чехия, Франция, Австрия) и Америку.

Влияние года на продуктивность

При поиске выбросов мы предполагали, что гетерогенность выборки достигается не только сортовыми различиями, но и тем, что растения в разные года различались климатические условия, что несомненно влияло на растения. Давайте проверим это гипотезу.

H0 - продуктивность не менялась от года к году

H1 - продуктивность отличалась в разные года

Сначала визуально оценим распределение

По построенным боксплотам видно, что среднее групп различается, поэтому гипотезу стоит проверить

model1 <- aov(productivity ~ year, data_without_na)
summary(model1)

##               Df  Sum Sq Mean Sq F value Pr(>F)    
## year           4  351522   87880   24.06 <2e-16 ***
## Residuals   1190 4346234    3652                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

По полученным графикам построенной модели, определим соблюдение допущений:

Полученные графики кажутся приемлемыми для использования модели и проведения post hoc теста.

В данном случае решили проверить в каких граппах есть значимые отличия с помощью SNK теста:

plot(agricolae::SNK.test(model1, trt = c('year')), main = "Результаты теста SNK по годам")

mtext("Продуктивность", side = 2, line = 3)  # Подпись оси Y
mtext("Год", side = 1, line = 2)

Наиболее сильно климатические условия на урожайность повлияли 2019 и в 2021, также образовалась отдельная группа, внутри которой статистически значимых влияний на продуктивность не обнаружено, это: 2018, 2017, 2020 года.

Построение моделей

Линейная модель

Возьмём все числовые переменные в качестве предикторов для предсказания продуктивности

model_prod <- lm(productivity ~ flowering_group+maturation_group+vegetation_period+protein_content+oil_content, data_without_na)
summary(model_prod)

## 
## Call:
## lm(formula = productivity ~ flowering_group + maturation_group + 
##     vegetation_period + protein_content + oil_content, data = data_without_na)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -137.041  -31.691   -0.026   31.938  177.910 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -135.2531    46.7716  -2.892 0.003901 ** 
## flowering_group.L   -49.7753    13.5435  -3.675 0.000248 ***
## flowering_group.Q    -1.8186    10.1731  -0.179 0.858152    
## flowering_group.C    -6.8416     8.1239  -0.842 0.399873    
## flowering_group^4   -11.5272     8.1703  -1.411 0.158546    
## flowering_group^5     5.5259    10.4461   0.529 0.596913    
## flowering_group^6    -6.3540     6.3646  -0.998 0.318327    
## flowering_group^7     0.0156     9.9404   0.002 0.998748    
## flowering_group^8     8.6407     9.3895   0.920 0.357631    
## maturation_group.L  -33.0893     9.3378  -3.544 0.000410 ***
## maturation_group.Q  -45.1036     6.7629  -6.669 3.94e-11 ***
## maturation_group.C    6.4134     4.4651   1.436 0.151172    
## maturation_group^4   -6.7338     2.7922  -2.412 0.016031 *  
## vegetation_period     1.9783     0.1654  11.958  < 2e-16 ***
## protein_content      -4.2524     0.6458  -6.584 6.86e-11 ***
## oil_content          12.6826     1.0111  12.544  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 45.86 on 1179 degrees of freedom
## Multiple R-squared:  0.4721, Adjusted R-squared:  0.4653 
## F-statistic: 70.28 on 15 and 1179 DF,  p-value: < 2.2e-16

Согласно значениям Pr(>|t|) все используемые предикторы являются важными. Но в любом случае попробуем проанализировать как измениться предсказательность модели при исключении одного/нескольких предикторов при обнаружении мультиколлинеарности между ними.

Посмотрим на графиках

par(mfrow=c(2,2))
plot(model_prod, which=1:4)

Оценим нормальность распределения остатков тестом Шапиро-Уилка

shapiro.test(residuals(model_prod))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_prod)
## W = 0.99869, p-value = 0.53

p-value = 0.1262 - остатки подчиняются нормальному распределению

Оценим коллинеарность предикторов с помощью corr.test

cor(select(data, vegetation_period, protein_content, oil_content))

##                   vegetation_period protein_content oil_content
## vegetation_period                 1              NA          NA
## protein_content                  NA               1          NA
## oil_content                      NA              NA           1

Многие из предикторов сильно коррелируют, например группа цветения c вегетационным периодом и группой созревания.

Оценим с помощью коэффициента инфляции дисперсии (VIF) из пакета car мультиколлениарность предикторов:

vif(model_prod)

##                        GVIF Df GVIF^(1/(2*Df))
## flowering_group    8.397932  8        1.142249
## maturation_group  10.424573  4        1.340471
## vegetation_period  2.223400  1        1.491107
## protein_content    1.618014  1        1.272012
## oil_content        1.718886  1        1.311063

Получили значения, не превышающие 5, однако попробуем убрать из ряда предикторов с наибольшим значением. Например, такие, как maturation_group и vegetation_period. Начнем с maturation_group

model_prod_1 <- update(model_prod, .~. - maturation_group) 
vif(model_prod_1)

##                       GVIF Df GVIF^(1/(2*Df))
## flowering_group   1.912530  8        1.041359
## vegetation_period 1.669979  1        1.292277
## protein_content   1.560378  1        1.249151
## oil_content       1.677473  1        1.295173

Значения коэффициентов vif немного снизились для остальных предикторов.

Попробуем оценить значимость всех предикторов с помощью drop1

drop1(model_prod, test = "F") 

## Single term deletions
## 
## Model:
## productivity ~ flowering_group + maturation_group + vegetation_period + 
##     protein_content + oil_content
##                   Df Sum of Sq     RSS    AIC  F value    Pr(>F)    
## <none>                         2480123 9159.3                       
## flowering_group    8     64676 2544799 9174.1   3.8432 0.0001775 ***
## maturation_group   4    202289 2682412 9245.0  24.0410 < 2.2e-16 ***
## vegetation_period  1    300804 2780926 9294.1 142.9959 < 2.2e-16 ***
## protein_content    1     91197 2571320 9200.5  43.3531 6.865e-11 ***
## oil_content        1    330983 2811106 9307.0 157.3424 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Согласно оценке, наибольшую значимость имеют предикторы: oil_content, vegetation_period. Наименее значимый - группа цветения (flowering_group).

Сравним модели с помощью информационного критерия Акаике (AIC), Adjusted R-squared и выберем наиболее подходящую.

1. Модель с исключением из предикторов группы созревания (maturation_group)

summary(model_prod_1)

## 
## Call:
## lm(formula = productivity ~ flowering_group + vegetation_period + 
##     protein_content + oil_content, data = data_without_na)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -137.558  -32.904    1.545   32.888  162.343 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -76.9648    46.3751  -1.660   0.0973 .  
## flowering_group.L -57.5438     9.8652  -5.833 7.02e-09 ***
## flowering_group.Q -37.5436     8.9469  -4.196 2.92e-05 ***
## flowering_group.C  -4.9303     8.3161  -0.593   0.5534    
## flowering_group^4  -1.0774     8.1919  -0.132   0.8954    
## flowering_group^5  -2.7247    10.5309  -0.259   0.7959    
## flowering_group^6   3.9505     6.0728   0.651   0.5155    
## flowering_group^7  -0.2820     9.8257  -0.029   0.9771    
## flowering_group^8  16.5387     9.4797   1.745   0.0813 .  
## vegetation_period   1.4209     0.1489   9.545  < 2e-16 ***
## protein_content    -4.6166     0.6585  -7.011 3.97e-12 ***
## oil_content        13.3814     1.0370  12.904  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 47.62 on 1183 degrees of freedom
## Multiple R-squared:  0.429,  Adjusted R-squared:  0.4237 
## F-statistic:  80.8 on 11 and 1183 DF,  p-value: < 2.2e-16

drop1(model_prod_1)

## Single term deletions
## 
## Model:
## productivity ~ flowering_group + vegetation_period + protein_content + 
##     oil_content
##                   Df Sum of Sq     RSS    AIC
## <none>                         2682412 9245.0
## flowering_group    8    367004 3049416 9382.2
## vegetation_period  1    206602 2889014 9331.7
## protein_content    1    111457 2793869 9291.7
## oil_content        1    377556 3059968 9400.4

Для модели 1 имеем:

Adjusted R-squared:  0.3918

AIC: 9302.5

2. Модель с исключением из предикторов группы цветения (flowering_group)

   model_prod_2 <- update(model_prod, .~. - flowering_group) 
   summary(model_prod_2)

   ## 
   ## Call:
   ## lm(formula = productivity ~ maturation_group + vegetation_period + 
   ##     protein_content + oil_content, data = data_without_na)
   ## 
   ## Residuals:
   ##      Min       1Q   Median       3Q      Max 
   ## -136.789  -32.833    0.273   31.376  166.425 
   ## 
   ## Coefficients:
   ##                     Estimate Std. Error t value Pr(>|t|)    
   ## (Intercept)        -133.0844    46.4649  -2.864  0.00425 ** 
   ## maturation_group.L  -57.3939     6.8717  -8.352  < 2e-16 ***
   ## maturation_group.Q  -46.9199     4.5602 -10.289  < 2e-16 ***
   ## maturation_group.C    0.3291     3.5852   0.092  0.92687    
   ## maturation_group^4   -4.7701     2.6328  -1.812  0.07027 .  
   ## vegetation_period     1.8183     0.1628  11.169  < 2e-16 ***
   ## protein_content      -4.1793     0.6484  -6.446 1.67e-10 ***
   ## oil_content          13.2361     0.9973  13.271  < 2e-16 ***
   ## ---
   ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   ## 
   ## Residual standard error: 46.3 on 1187 degrees of freedom
   ## Multiple R-squared:  0.4583, Adjusted R-squared:  0.4551 
   ## F-statistic: 143.5 on 7 and 1187 DF,  p-value: < 2.2e-16

   drop1(model_prod_2)

   ## Single term deletions
   ## 
   ## Model:
   ## productivity ~ maturation_group + vegetation_period + protein_content + 
   ##     oil_content
   ##                   Df Sum of Sq     RSS    AIC
   ## <none>                         2544799 9174.1
   ## maturation_group   4    504616 3049416 9382.2
   ## vegetation_period  1    267435 2812234 9291.5
   ## protein_content    1     89077 2633876 9213.2
   ## oil_content        1    377609 2922408 9337.4

Для модели 2 имеем:
Adjusted R-squared:  0.4033
 AIC: 9279.7

3. Модель с исключением из предикторов данных о вегетационном периоде (vegetation_period)

   model_prod_3 <- update(model_prod, .~. -vegetation_period) 
   summary(model_prod_3)

   ## 
   ## Call:
   ## lm(formula = productivity ~ flowering_group + maturation_group + 
   ##     protein_content + oil_content, data = data_without_na)
   ## 
   ## Residuals:
   ##      Min       1Q   Median       3Q      Max 
   ## -153.765  -32.527    1.365   34.225  164.032 
   ## 
   ## Coefficients:
   ##                    Estimate Std. Error t value Pr(>|t|)    
   ## (Intercept)        100.4945    44.8931   2.239   0.0254 *  
   ## flowering_group.L  -33.0707    14.2588  -2.319   0.0205 *  
   ## flowering_group.Q   -2.5405    10.7676  -0.236   0.8135    
   ## flowering_group.C   -5.7597     8.5983  -0.670   0.5031    
   ## flowering_group^4  -12.7773     8.6472  -1.478   0.1398    
   ## flowering_group^5    7.3222    11.0556   0.662   0.5079    
   ## flowering_group^6   -8.4767     6.7341  -1.259   0.2084    
   ## flowering_group^7    7.3414    10.5015   0.699   0.4846    
   ## flowering_group^8    2.5783     9.9239   0.260   0.7951    
   ## maturation_group.L  13.7320     8.9728   1.530   0.1262    
   ## maturation_group.Q -46.7610     7.1567  -6.534 9.51e-11 ***
   ## maturation_group.C   9.6968     4.7172   2.056   0.0400 *  
   ## maturation_group^4  -6.3021     2.9551  -2.133   0.0332 *  
   ## protein_content     -5.1364     0.6791  -7.563 7.87e-14 ***
   ## oil_content         12.9425     1.0699  12.096  < 2e-16 ***
   ## ---
   ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   ## 
   ## Residual standard error: 48.55 on 1180 degrees of freedom
   ## Multiple R-squared:  0.408,  Adjusted R-squared:  0.401 
   ## F-statistic:  58.1 on 14 and 1180 DF,  p-value: < 2.2e-16

   drop1(model_prod_3)

   ## Single term deletions
   ## 
   ## Model:
   ## productivity ~ flowering_group + maturation_group + protein_content + 
   ##     oil_content
   ##                  Df Sum of Sq     RSS    AIC
   ## <none>                        2780926 9294.1
   ## flowering_group   8     31308 2812234 9291.5
   ## maturation_group  4    108088 2889014 9331.7
   ## protein_content   1    134817 2915743 9348.7
   ## oil_content       1    344846 3125772 9431.8

Для модели 3 имеем:

Adjusted R-squared:  0.3512
AIC: 9379.7

При этом для исходной модели с сохранением всех предиктров:

Adjusted R-squared:  0.4141
AIC: 9258.8

Таким образом, точность предсказания исходной модели максимальная в сравнении с тремя рассмотренными выше и составляет 41.41%; также для нее информационный критерий Акаике был минимален

Для предсказания продуктивности сортов наиболее подходящей будет модель, учитывающая группу цветения, вегетационный период, группу созревания, содержание белков и масел

Смешанная линейная модель

model5 <- lme4::lmer(productivity ~ leaf_shape + maturation_group + lodging_type + (1 | year), data = data_without_na)
summary(model5)

## Linear mixed model fit by REML ['lmerMod']
## Formula: productivity ~ leaf_shape + maturation_group + lodging_type +  
##     (1 | year)
##    Data: data_without_na
## 
## REML criterion at convergence: 12927.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.00933 -0.69312 -0.01955  0.68005  2.88914 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  year     (Intercept)  476.5   21.83   
##  Residual             3004.9   54.82   
## Number of obs: 1195, groups:  year, 5
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)         151.888     12.917  11.759
## leaf_shaperound       6.679      5.703   1.171
## maturation_group.L   -1.175      6.025  -0.195
## maturation_group.Q  -71.319      5.231 -13.635
## maturation_group.C   12.295      4.244   2.897
## maturation_group^4   -6.904      3.114  -2.217
## lodging_typeno       -4.654      6.832  -0.681
## lodging_typeyes     -50.538     12.476  -4.051
## 
## Correlation of Fixed Effects:
##             (Intr) lf_shp mtr_.L mtr_.Q mtr_.C mtr_^4 ldgng_typn
## leaf_shprnd -0.381                                              
## mtrtn_grp.L -0.059  0.001                                       
## mtrtn_grp.Q  0.155 -0.069 -0.299                                
## mtrtn_grp.C -0.105  0.129  0.477 -0.221                         
## mtrtn_grp^4  0.077 -0.092 -0.116  0.381 -0.090                  
## lodgng_typn -0.481 -0.054  0.035 -0.041  0.047 -0.042           
## ldgng_typys -0.279  0.002 -0.097 -0.103 -0.045 -0.057  0.512

## [[1]]

## `geom_smooth()` using formula = 'y ~ x'

## 
## [[2]]
## [[2]]$year

## `geom_smooth()` using formula = 'y ~ x'

## 
## 
## [[3]]

## 
## [[4]]

## `geom_smooth()` using formula = 'y ~ x'

• Scaled residuals показывают распределение остатков модели. Нормальное распределение остатков позволяет предположить, что модель адекватна.

  Случайные эффекты:

• year имеет вариацию 464.3 и стандартное отклонение 21.55, что говорит о значительном влиянии года на продуктивность.

  Фиксированные эффекты:

• Intercept: Оценка пересечения равна 202.435. Это базовый уровень продуктивности для группы, когда все остальные переменные равны нулю.

• leaf_shaperound: Небольшой положительный эффект (+2.670), но он не статистически значим (t value = 0.439).

• maturation_group: Отрицательный значимый эффект (-8.021), что указывает на то, что увеличение группы созреваемости связано с уменьшением продуктивности.

• lodging_typeno: Отрицательный эффект (-7.796) при отсутствии уклона, но не статистически значимый.

• lodging_typeyes: Заметный отрицательный эффект (-67.445) с высоким t-значением (-5.053), что указывает на значительное снижение продуктивности, когда растения предрасположены к уклону.

Выводы из модели показывают, что maturation_group и lodging_type имеют значительное влияние на продуктивность. Возможная предрасположенность к уклону существенно снижает продуктивность.

Вывод

1. С помощью дисперсионного анализа ANOVA мы установили, что регион происхождения влияет на продуктивность: отличия между местами произрастания замечены в странах находящихся в разных частях света и соотвественно различных климатических зонах.

2. С помощью дисперсионного анализа ANOVA мы установили, что продуктивность сои также отличается в зависимости от года, с чем связана часть гетерогенности и выбросов в нашей выборке. Наиболее сильно климатические условия на урожайность повлияли 2019 и в 2021 году.

3. Подобрали модель предсказывающую продуктивность которая в качестве предикторов учитывает группу цветения, вегетационный период, группу созревания, содержание белков и масел.